*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0

* This do file generates internal and external distance from prime location boundaries in US MSAs
* These distances will be used in the BDD plots that show how emplyoment densities change across the prime location boundary
 
* Generate temporary working directory
	qui capture mkdir "$temp/TEMPDISTBORDER" // Generate temporary working directory
 
* Start loop over pvalues *****************************************************
	foreach x in $plist  { 
			   
	* set up metro list
		qui u "$data_USMETROS/Raw Numeric Data/METRO LIST/METROS.dta", clear	
		display "...Start working on pvalue `x'..."
		qui tab cbsafp
		qui local Ncbsafp = r(N) // Total number of cities in list
		qui levelsof  cbsafp, local(MID)	
		local count=1
	
	*  Loop over MSAs starts ***************************************************
		foreach m_ID of local MID {
			display "...working on pval `x' and MSA `m_ID'; `count' of `Ncbsafp'..."
		* Load grid data
			qui use "$temp/grid_PL_output_p`x'/PL_grids_p`x'_`m_ID'", clear
			preserve
			qui keep square_id lat lon cell_id
			qui rename lat lat_p
			qui rename lon lon_p
			qui	rename square_id square_id_p														
			qui save "$temp/TEMPDISTBORDER/points" ,replace 
			restore 
		* compute distances						
			preserve 
			qui geonear square_id lat lon  using  "$temp/TEMPDISTBORDER/points", n(square_id_p lat_p lon_p) ignoreself within(.354) long
			qui collapse (count) km_ , by(square_id)
			qui rename km_ neighbors
			qui gen border_cell=1 if neighbors!=8
			qui gen border_cell_corner=1 if neighbors==7
			qui replace border_cell=0 if border_cell==.
			qui keep square_id border_cell border_cell_corner 
			qui	save "$temp/TEMPDISTBORDER/neighbors_ids" ,replace 
			restore 
  		* Merge IDs			   
			qui merge 1:1 square_id using "$temp/TEMPDISTBORDER/neighbors_ids", nogen
		* Calculate internal distances 
			preserve 
			qui keep if border_cell==1
			qui keep square_id lat lon cell_id
			qui rename square_id nborid
			qui rename lat nborlat
			qui rename lon nborlon
			qui save "$temp/TEMPDISTBORDER/nbordata_border", replace  
			restore 
			qui	geonear  square_id lat lon using "$temp/TEMPDISTBORDER/nbordata_border" , ///
				genstub(borderdist_) n(nborid nborlat nborlon) near(1) 
		* add internal distance from centroid because above is measured centroid-centroid 
		* but we only want centroid-border
			qui gen border_dist_internal=-(km_to_borderdist_+.125)
		*replace for those that are only corner 
			qui replace border_dist_internal=-(sqrt(.125^2+.125^2)) if  border_cell_corner==1
		* clear up variables 
			qui drop borderdist_ km_to_borderdist_   
			qui keep square_id   border_dist_internal
			qui save "$temp/TEMPDISTBORDER/internal_distances", replace
		*** calculate external distances
		* join full city grid
			qui use "$temp/grid_PL_output_p`x'/grid_PL_output_p`x'_`m_ID'", clear
			qui keep  lat lon  square_id PL cell_id
			qui geonear  square_id lat lon using "$temp/TEMPDISTBORDER/nbordata_border"  , genstub(BD_external) n(nborid nborlat nborlon) near(1)  
			qui sort km_to_BD_external
			qui merge 1:1 square_id using "$temp/TEMPDISTBORDER/internal_distances", nogen
		*substract   125 meters as we are calculating centroid to centroid distances here
			qui gen border_dist_external=km_to_BD_external-(.125) if PL==0 
			qui gen distance_toPLborder=border_dist_internal if PL==1
			qui replace distance_toPLborder=border_dist_external if PL==0
			qui keep   lat lon cell_id square_id distance_toPLborder
			qui rename distance_toPLborder dist_`x'
 			qui save "$temp/TEMPDISTBORDER/DIST_`m_ID'_`x'", replace 
		* Update MSA counter
			local count = `count'+1
		}
		* Loop over MSA ended **************************************************
	
	* Clear memoray and append all MSAs	
		display "...finalizing pvalue `x'..."
		clear
		* Loop over MSA starts 
		qui foreach m_ID of local MID {
			append using "$temp/TEMPDISTBORDER/DIST_`m_ID'_`x'"
			save "$temp/TEMPDISTBORDER/ALLDIST_`x'", replace 
		}
	}
		* Loop over MSAs ends
* Loop over pvalues ended ******************************************************

* Finalize distances data set 
	display "..finalizing distances data set..."
	clear
 * Merge with whole grid data and export
	qui use "$temp/grid_PL_output_p98_5/grid_PL_output_p98_5__all.dta", clear 
	qui keep cell_MFG_emp cell_PS_emp cell_TS_emp cell_O_emp cell_ST_emp emp area_g developable cbsafp cell_id square_id
	qui foreach x in  $plist {  
		merge 1:1 cell_id  using "$temp/TEMPDISTBORDER/ALLDIST_`x'",  keepusing(dist_`x') nogen    update 
		label var dist_`x' "distance to border of PL (internal distances negative) "
	 }
	*label var distance_toPLborder "distance to border of PL (internal distances negative)"
		save  "$temp/USMETROS_PL_distances", replace

* Script ends	
